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HIGHLIGHTS 


►  An  LTI  model,  in  the  form  of  a  state  space  model,  is  proposed  for  solid-phase  diffusion  in  physics-based  lithium  ion  cell  models. 

►  The  proposed  model  can  be  used  for  spherical  and  non-spherical  particles. 

►  Impact  of  different  particle  shapes  on  electrochemical  performances  can  be  investigated. 

►  The  model  is  much  faster  than  solving  the  full  model. 
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Physics-based  lithium  ion  models  are  widely  used  to  predict  the  electrochemical  behavior  of  lithium  ion 
cells.  The  implementation  of  such  a  model  typically  requires  solving  a  diffusion  problem  in  solid  parti¬ 
cles.  A  linear  time-invariant  (LTI)  model  is  proposed  for  the  solid-phase  diffusion  problem.  This  LTI  model 
can  be  used  for  spherical  and  non-spherical  particles.  For  spherical  particles,  results  from  using  the  LTI 
model  are  compared  with  those  from  solving  full  diffusion  equation,  and  excellent  agreement  is  ach¬ 
ieved.  The  LTI  model  solves  only  a  few  equations,  and  thus  it  runs  much  faster  than  the  model  solving  the 
full  diffusion  equation.  Impact  of  particle  shapes  on  the  electrochemical  behavior  is  investigated  after  the 
model  is  validated. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

The  lithium  ion  battery  is  a  preferred  candidate  as  a  power 
source  for  hybrid  electric  vehicle  (HEV)  and  electric  vehicle  (EV) 
due  to  its  outstanding  characteristics  such  as  high  energy 
density,  high  voltage,  low  self-discharge  rate,  and  good  stability 
among  others.  Physics-based  lithium  ion  models  are  widely  used 
to  predict  the  electrochemical  behavior  of  lithium  ion  cells  [1—3], 
The  implementation  of  such  a  model  typically  requires  solving 
a  diffusion  problem  in  solid  particles,  which  are  used  to  model 
the  porous  electrodes.  Due  to  the  large  number  of  particles 
involved,  the  particle  shapes  are  assumed  to  be  spherical  so  that 
the  diffusion  equation  becomes  one  dimensional.  However,  it  is 
still  computationally  expensive  to  solve  a  large  number  of 
Id  diffusion  equations.  For  particles  of  arbitrary  shapes  or 
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more  general  porous  structures,  three  dimensional  diffusion 
equations  are  required  to  characterize  the  diffusion  process 
in  such  structures.  A  physics-based  model,  which  involves 
solving  3d  diffusion  equations,  becomes  prohibitively  expensive 
computationally. 

For  spherical  particles,  several  methods  have  been  proposed  to 
reduce  the  size  of  the  problem.  They  belong  to  either  the  global 
approach  or  the  local  approach.  In  the  global  approach,  the  entire 
model  is  reduced.  Such  a  global  approach  is  used  in  [4]  through 
proper  orthogonal  decomposition  (POD).  In  such  an  approach,  the 
space  of  the  solution  is  estimated  by  testing  system  responses 
from  typical  boundary  conditions  to  a  full  model.  A  set  of  basis  are 
then  formed  for  the  solution  space  in  such  a  way  that  the  leading 
coordinates  have  much  larger  magnitude  and  thus  the  rest  of  the 
coordinates  can  be  truncated.  The  reduced  model  is  then  formu¬ 
lated  based  on  the  new  truncated  basis.  Because  the  physics-based 
model  is  highly  non-linear,  when  actual  boundary  conditions 
differ  from  the  typical  boundary  conditions  used  to  identify  the 
reduced  model  the  error  could  be  large.  This  approach  cannot  be 
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extended  to  non-spherical  particles  easily  since  a  full  model  is 
needed  in  the  first  place  before  a  reduced  model  could  be  gener¬ 
ated.  Another  global  approach  is  proposed  in  [5].  In  this  approach, 
a  full  model  is  linearized  and  transfer  functions  are  created  for  the 
linearized  model.  Errors  could  be  large  in  using  such  an  approach 
if  perturbation  is  large  since  linearization  assumes  small  pertur¬ 
bation.  The  method  also  relies  on  spherical  particles  for  an 
analytical  transfer  function  and  thus  cannot  be  extended  to  non- 
spherical  particles. 

Another  popular  approach  is  the  local  approach,  in  which  only 
the  linear  diffusion  equation  is  reduced  and  the  non-linear  part  of 
the  model  is  intact.  Different  local  methods  differ  in  how  to  obtain 
approximate  solutions  to  the  diffusion  problem  for  a  spherical 
particle.  Approximation  can  be  performed  either  in  the  time- 
domain  or  in  the  frequency-domain.  One  time-domain  approach 
is  to  take  advantage  of  superposition  [1,2],  Another  time-domain 
approach  makes  the  assumption  that  the  concentration  within 
each  spherical  particle  can  be  approximated  with  a  parabolic 
profile  [6,7],  A  relatively  recent  time-domain  approach  is  to  obtain 
an  approximate  solution  by  truncating  the  analytical  solution  of  an 
infinite  series  and  then  adding  an  estimate  term  for  the  truncation 
error  [8],  In  all  of  the  time-domain  local  methods  mentioned, 
spherical  particles  are  assumed.  In  frequency-domain  local 
approximation  methods,  a  state  space  model  is  created  to 
approximate  the  transfer  function  of  the  system.  In  Ref.  [9],  high- 
order  poles  of  the  transfer  function  are  truncated  and  the  lower- 
order  poles  are  grouped  together  and  approximated  using  a  state 
space  model.  In  Ref.  [10],  a  discrete-time  state  space  model  is 
derived  from  a  known  transfer  function.  In  Refs.  [9,10],  analytical 
transfer  functions  are  assumed  to  be  known,  and  thus  both 
methods  are  limited  to  spherical  particles. 

In  this  paper,  an  LTI  model,  a  frequency-domain  local  approach, 
is  proposed  to  accurately  model  the  diffusion  process  for  spherical 
and  non-spherical  particles.  This  method  shares  similarity  with  [9] 
and  [10]  in  that  it  also  seeks  a  state  space  model  that  approximates 
the  transfer  function  of  the  system.  However,  the  current  method 
does  not  assume  a  known  transfer  function.  The  transfer  function  is 
calculated  numerically  and  thus  it  can  be  applied  to  non-spherical 
particles.  The  way  to  identify  the  state  space  model  is  also  different. 
The  method  used  in  Ref.  [9]  does  not  give  the  same  accuracy  as  the 
LTI  method  in  this  paper.  Compared  with  Ref.  [10],  the  current 
approach  obtains  the  continuous-time  state  space  model  directly 
rather  than  having  to  obtain  the  discrete-time  state  space  model 
first.  The  LTI  method  has  been  widely  used  to  model  devices  and 
subsystems  for  the  purpose  of  transient  analysis  in  power  systems 
[11,12],  signal  integrity  characterization  of  microwave  systems 
[13,14],  and  battery  thermal  systems  [15,16],  In  applying  the  LTI 
method,  one  essentially  obtains  the  transfer  function  of  the  system 
numerically  (or  by  testing  in  some  cases)  and  then  a  state  space 
model  is  identified  to  have  the  same  transfer  function  using  the 
vector  fitting  method  [17],  Afterward,  the  state  space  model  can  be 
used  to  simulate  the  system  in  the  time-domain.  The  LTI  approach 
can  give  highly  accurate  results  and  yet  the  size  of  the  model  can  be 
very  small  if  the  state  space  model  is  properly  identified.  Note  that 
the  LTI  approach  also  relies  on  linearity,  but  it  does  not  require  the 
existence  of  an  analytical  solution  since  the  transfer  function  of  the 
system  is  obtained  numerically  (or  by  testing).  For  the  current 
application,  this  implies  that  the  method  can  be  used  for  particles 
of  any  shape. 

The  paper  is  organized  as  follows.  Section  2  describes  the  LTI 
approach  for  the  solid-phase  diffusion  problem  in  particles.  Section 
3  integrates  the  LTI  model  into  the  physics-based  model  using 
spherical  particles  and  validates  its  results  against  a  full  model 
solving  particle  diffusion  equations  directly.  In  Section  4,  non- 
spherical  particles  are  used  and  the  impact  of  different  shapes  of 


particles  on  the  electrochemical  behavior  is  investigated.  Finally, 
Section  5  is  the  conclusion. 

2.  LTI  modeling  for  particle  solid-phase  diffusion 

In  using  the  LTI  approach,  the  problem  is  treated  like  a  system.  A 
system  is  an  entity  that  processes  a  set  of  input  signals  (or  simply 
called  inputs)  and  yields  another  set  of  output  signals  (or  simply 
called  outputs).  In  such  a  system  view,  only  the  input/output  rela¬ 
tionship  of  the  system  is  of  interest  to  the  user  and  the  inner 
structure  of  the  system  is  not.  In  the  solid-phase  diffusion  problem, 
the  molar  flux  at  the  particle  surface  as  a  function  of  time  is  used  as 
the  input  and  the  surface-averaged  concentration  increase  at  the 
particle  surface  as  a  function  of  time  is  used  as  the  output.  We  are 
interested  in  the  relationship  between  the  surface  molar  flux  and 
averaged  surface  concentration  increase.  Concentration  distribu¬ 
tion  inside  a  particle  is  not  of  great  interest.  Solid-phase  diffusion 
problem  so  described  is  a  system.  More  importantly,  it  is  not  only 
a  system,  but  it  is  also  an  LTI  system.  Any  state  space  model  is  also 
an  LTI  system.  One  important  feature  about  LTI  systems  is  that  if 
two  LTI  systems  have  the  same  impulse  or  step  response,  the  two 
systems  behave  identically  in  that  the  outputs  of  the  two  systems 
are  the  same  provided  that  the  inputs  to  the  two  systems  are  the 
same.  This  feature  allows  us  to  use  the  state  space  model  to 
simulate  the  solid-phase  diffusion  problem  provided  that  its  step 
response  is  curve-fitted  to  that  of  the  solid-phase  diffusion 
problem.  Note  that  if  two  LTI  systems  have  the  same  impulse  or 
step  response,  they  have  the  same  transfer  function.  So,  when 
matching  the  step  response  of  two  LTI  systems,  effectively  we  are 
matching  the  transfer  function  of  the  two  LTI  systems.  Since  the  size 
of  the  state  space  model  is  small,  it  runs  very  fast  compared  with 
the  full  model  which  solves  the  complete  diffusion  equation 
directly  using  numerical  methods. 

The  diffusion  of  lithium  ions  in  a  solid-phase  particle  follows  the 
Fick’s  law  and  is  described  by  the  following  partial  differential 
equations: 

ft  =  DV2c  (1) 


(2) 


-D^-rnm  (3) 

8nlo 

where  c  is  the  concentration  of  lithium  in  the  solid  particle;  Co 
is  the  initial  concentration  of  lithium  in  the  solid  particle,  assumed 
to  be  a  constant;  D  is  the  solid-phase  diffusion  coefficient  for 
lithium  in  the  particle,  assumed  to  be  constant;  j(t)  is  the 
boundary  molar  flux,  assumed  to  be  a  function  of  time  only  for 
each  particle;  and  Q  denotes  the  boundary  of  the  particle,  which  is 
assumed  to  fixed.  We  are  interested  in  the  relationship  between 
j(t)  and  averaged  surface  concentration  increase  of  the  system. 
Since  the  diffusion  equation  and  the  boundary  condition  are  linear 
and  diffusivity  is  a  constant,  the  system  so  described  is  an  LTI 
system. 

A  general  state  space  model  is  typically  written  as  follows: 


y  = 


Ax  +  Bu 
Cx  +  Du 


(4) 


x  is  the  state  vector;  y  is  the  output  vector;  u  is  the  input  vector;  A,  B, 
C,  D  are  constant  coefficient  matrices  of  proper  sizes.  Since  Eqn.  (4) 
is  linear  with  constant  coefficients,  the  system  is  an  LTI  system.  For 
the  state  space  model  used  to  simulate  the  solid-phase  diffusion 
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problem,  x  has  no  physical  meaning;  y  is  a  scalar  representing 
averaged  surface  concentration  increase;  u  is  a  scalar  representing 
the  transient  flux  j(t).  A  critical  step  in  LTI  modeling  is  to  identify 
matrices  A,  B,  C,  D  such  that  the  state  space  model  gives  the  same 
step  response  as  the  original  diffusion  problem  governed  by  Eqn. 
(1)— (3).  Once  that  is  accomplished,  characteristics  of  LTI  systems 
guarantee  that  the  state  space  model  can  be  used  to  replace  Eqn. 
(1 ) — (3)  to  model  the  relationship  between  the  surface  molar  flux 
and  averaged  surface  concentration  increase  with  excellent 
accuracy. 

To  demonstrate  the  LTI  approach,  state  space  models  are  con¬ 
structed  for  three  different  solid-phase  structures,  a  spherical 
particle,  an  elliptical  particle,  and  a  more  general  porous  structure. 
In  the  first  example,  a  spherical  particle  is  used.  The  first  step  in  LTI 
model  identification  is  to  obtain  the  step  response  of  the  system 
modeled.  In  order  to  obtain  the  step  response  of  the  particle 
diffusion  problem,  the  diffusion  problem  is  solved  numerically 
using  FLUENT,  a  CFD  code  from  ANSYS.  The  concentration  distri¬ 
bution  at  the  end  of  the  step  response  simulation  is  shown  in  Fig.  1 
along  with  the  geometry  and  mesh  used  in  the  CFD  model.  Note 
that  any  CFD  code  can  be  used  to  solve  the  particle  diffusion 
problem  by  using  analogy  between  thermal  diffusion  and  species 
diffusion. 

After  the  diffusion  problem  is  solved  for  its  step  response,  a  state 
space  model  of  Eqn.  (4)  is  identified  by  curve-fitting  its  step 


response  to  that  of  the  CFD  model.  The  curve-fitting  was  performed 
in  the  frequency-domain  using  the  vector  fitting  (VF)  method  [17], 
A  brief  introduction  of  the  VF  method  applied  to  the  diffusion 
problem  is  provided  in  Appendix  A.  The  curve-fitting  results  are 
shown  in  Fig.  2.  It  can  be  seen  that  the  state  space  model  gives  the 
same  step  response  as  the  CFD  model  solving  Eqn.  (1)— (3).  The 
excellent  accuracy  in  Fig.  2  indicates  that  the  curve-fitting  by  the  VF 
method  is  very  accurate  for  this  diffusion  problem.  The  state  space 
model  has  an  order  of  7,  and  so  only  7  ordinary  differential  equa¬ 
tions  are  solved  in  the  state  space  model.  Higher  order  state  space 
models  can  be  used  for  greater  accuracy.  But  from  Fig.  2,  it  shows 
that  an  order  of  7  is  sufficient,  and  validation  in  Section  3  also 
indicates  that  7th  order  is  accurate  enough  when  such  a  model  is 
integrated  into  physics-based  lithium-ion  cell  models.  After  the 
state  space  model  is  identified,  it  can  be  used  to  simulate  the 
diffusion  problem  under  any  transient  flux  boundary  rather  than 
just  step  flux  boundary. 

For  a  spherical  particle,  its  step  response  can  be  obtained 
analytically  [8],  So,  the  state  space  model  could  also  be  identified 
using  the  analytical  solution.  Such  fitting  results  and  its  corre¬ 
sponding  state  space  model  are  provided  in  Appendix  B  for 
reference. 

Note  that  the  model  generation  process  of  the  LTI  approach  does 
not  depend  on  the  shape  of  the  particles.  In  the  second  example,  an 
elliptical  particle  is  used.  The  model  generation  process  is  the  same 


Fig.  1.  Geometry  and  mesh  for  a  spherical  particle  and  its  concentration  distribution  at  the  end  of  a  step  response  run. 
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Fig.  2.  Step  responses  from  the  CFD  model  and  the  LTI  model  for  the  spherical  particle.  Fig.  4.  Step  responses  from  the  CFD  model  and  the  LTI  model  for  the  elliptical  particle. 


as  before.  The  concentration  distribution  at  the  end  of  the  step 
response  simulation  along  with  the  geometry  and  mesh  used  in  the 
CFD  model  is  shown  in  Fig.  3.  Step  responses  from  the  CFD  model 
and  the  LTI  model  are  compared  in  Fig.  4.  The  excellent  agreement 


shown  in  Fig.  4  indicates  that  curve-fitting  performed  by  the  VF 
method  is  excellent.  This  also  implies  that  subsequent  simulation 
by  the  LTI  model  under  any  transient  flux  j(t)  will  give  excellent 
accuracy. 


Fig.  3.  Geometry  and  mesh  for  an  elliptical  particle  and  its  concentration  distribution 


i  end  of  a  step  response  run. 


X  Hu  et  al.  /Journal  of  Power  Sources  214  (2012)  40-50 


The  above  elliptical  example  showed  that  the  model  generation 
process  of  the  LTI  approach  is  quite  general.  As  a  matter  of  fact,  the 
porous  structure  of  electrodes  does  not  even  have  to  be  modeled 
using  particles.  Instead,  any  porous  structure  shape  could  be  used 
without  simplification.  In  this  third  example,  a  porous  structure  is 
used,  which  consists  of  spherical  particles  packed  together  with 
overlapping.  For  simplicity,  this  porous  structure  will  be  referred  to 
as  a  porous  particle  in  the  rest  of  the  paper.  For  this  porous  particle, 
the  LTI  model  identification  process  is  the  same  as  before.  The 
concentration  distribution  at  the  end  of  the  step  response  simula¬ 
tion  along  with  the  geometry  and  mesh  used  in  the  CFD  model  is 
shown  in  Fig.  5.  Step  responses  from  the  CFD  model  and  the  LTI 


model  are  compared  in  Fig.  6.  Since  concentration  changes  rapidly 
close  to  time  of  zero,  a  plot  using  log  time  scale  is  also  shown  in  Fig.  6. 
Both  plots  in  Fig.  6  showed  excellent  accuracy.  High  accuracy  near 
time  of  zero  is  important  for  accurate  prediction  of  transient 
behavior  of  a  battery  cell  as  shown  shortly.  Log  time  scale  plots  for 
the  spherical  and  elliptical  particles  are  also  excellent.  For  simplicity, 
only  the  log  time  scale  plot  for  the  porous  particle  is  shown. 

The  claim,  that  an  LTI  model  could  be  used  to  simulate  the 
diffusion  problem  under  any  transient  flux  of  j(t)  provided  that  its 
step  response  matches  well  with  that  of  the  CFD  model,  is  verified 
next.  To  do  that,  a  rather  arbitrary  flux  of  j(t)  shown  in  Fig.  7  is  used 
as  the  input  to  the  CFD  model  and  the  corresponding  LTI  model  for 


X  Hu  et  al.  /  Journal  of  Power  Sources  214  (2012)  40-50 


15030 
_  15010 
E  14990 
E.  14970 
.2  14950 
1  14930 


14870 


100  200  300  400 

Time  (sec) 


14850 


Time  (sec) 


Fig.  6.  Step  responses  from  the  CFD  model  and  the  LTI  model  for  the  porous  particle,  a) 
Regular  scale  on  time,  b)  Log  scale  on  time. 


the  porous  particle.  And  their  surface  concentration  results  are  then 
compared  in  Fig.  8.  It  is  clear  from  Fig.  8  that  the  state  space  model 
gives  excellent  results  under  such  an  arbitrary  transient  flux  of  j{t). 
Therefore,  the  claim  is  verified. 

If  a  spherical  particle  and  an  elliptical  particle  are  created  to 
have  the  same  volume,  the  difference  in  surface  area  would  cause 
different  surface  concentration  behavior.  The  elliptical  particle  has 
more  surface  area  and  therefore  under  the  same  surface  flux  its 
surface  concentration  increases  more  quickly.  On  the  other  hand,  if 
the  same  total  flux,  namely  the  surface  integration  of  j(t),  is  applied 


Time  (sec) 

Fig.  8.  Surface  concentration  from  the  CFD  model  and  the  LTI  model  under  j(t)  shown 
in  Fig.  7  for  the  porous  particle. 


to  the  two  particles,  the  elliptical  particle  experiences  smaller  flux 
and  thus  its  surface  concentration  increases  less  rapidly.  The  second 
scenario  can  be  used  to  discuss  the  impact  of  different  particle 
shapes  on  battery  performance.  The  same  argument  above  also 
applies  to  particles  of  different  volumes  and  porous  particles  after 
proper  scaling.  The  porous  particle  used  in  this  paper,  which  is 
generated  using  spherical  particles  of  the  same  size  but  with 
overlapping,  has  the  opposite  effect  compared  with  the  elliptical 
particle.  This  is  because  the  porous  particle  has  less  surface  area 
compared  with  the  spherical  particle  per  unit  volume.  Fig.  9  shows 
surface  concentration  of  the  three  different  particles  under  the 
same  total  flux,  and  the  different  behavior  is  what  was  just  dis¬ 
cussed.  Larger  surface  area  from  the  elliptical  particle  would  make 
such  a  particle  less  prone  to  saturation  and  depletion.  On  the  other 
hand,  smaller  surface  area  from  the  porous  particle  would  make 
such  a  particle  more  prone  to  saturation  and  depletion.  Impact  of 
different  particle  shapes  on  battery  performance  is  discussed  in 
details  in  Section  4. 


3.  LTI  model  validation 

In  Section  2,  LTI  models  are  created  and  validated  in  isolated 
environments.  In  this  section,  the  LTI  model  is  validated  when 
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Fig.  10.  Cell  potential  vs.  capacity  at  different  discharge  rates  using  the  full  model 
(denoted  by  solid  lines)  and  the  reduced  model  (denoted  by  dotted  lines). 


integrated  into  the  physics-based  lithium  ion  cell  model.  Such 
a  model  will  be  referred  to  as  the  reduced  model.  For  comparison, 
a  model  solving  the  complete  solid-phase  diffusion  equation 
directly  is  created.  This  model  will  be  referred  to  as  the  full  model. 
Spherical  particles  are  used  since  the  full  model  can  only  allow  for 
spherical  particles.  Properties  and  dimensions  of  the  two  models 
are  otherwise  identical  and  they  are  all  from  Ref.  [4],  Both  models 
are  created  using  the  VHDL-AMS  modeling  language  and  solved  in 
Simplorer,  a  system  simulator  from  ANSYS.  Generation  of  physics- 
based  lithium  ion  cell  models  using  VHDL-AMS  modeling 
language  is  discussed  in  detail  in  Ref.  [18], 

Fig.  10  showed  the  cell  discharge  curves  using  the  full  model  and 
the  reduced  model.  It  can  be  seen  from  Fig.  10  that  excellent  results 
are  obtained  by  using  the  reduced  model  even  at  a  discharge  rate  of 
10  C.  For  10  C,  the  full  model  needs  200  equations  for  each  particle 
to  obtain  accurate  results  when  spatial  discretization  is  uniform, 
and  for  the  rest  of  the  discharge  curves,  100  equations  per  particle  is 
used  to  obtain  adequate  accuracy  using  uniform  spatial  dis¬ 
cretization.  Fig.  11  showed  the  comparison  of  concentration  at  four 
specified  interfaces.  Again,  excellent  results  are  achieved  by  using 
the  reduced  model.  Fig.  12  compares  the  flux  j(t)  at  four  particle 


Time  (sec) 


Time  (sec) 

Fig.  12.  Surface  molar  flux  j(t)  obtained  from  the  full  model  (denoted  by  solid  lines) 
and  those  from  the  reduced  model  (denoted  by  dotted  lines)  at  1  C  discharge  rate  for 
a  few  selected  particles. 


surfaces,  and  Fig.  13  compares  the  surface  concentration  at  the 
same  four  particles.  Excellent  results  are  achieved  by  using  the 
reduced  model  in  both  Figs.  12  and  13.  Note  that  at  the  start  of  the 
discharging  surface  concentration  and  flux  change  rapidly  because 
of  sudden  change  of  current  from  0  to  1  C  at  time  of  zero.  And  the 
reduced  model  can  capture  all  the  transient  behavior  at  the  particle 
surface.  In  Fig.  14,  a  charge  discharge  cycle  is  simulated  for 
comparison.  Again,  we  observe  excellent  results  from  the  reduced 
model. 

All  of  the  comparisons  performed  above  show  that  the  reduced 
model  can  accurately  model  the  cell  behavior  and  thus  can  be  used 
to  replace  the  full  model.  While  the  full  model  solves  for  100 
equations  for  each  particle  (except  for  200  being  used  for  10  C 
discharge),  the  reduced  model  solves  only  7  equations  per  particle. 
Because  of  that,  the  reduced  model  runs  much  faster.  For  the  cell¬ 
cycling  simulation  shown  in  Fig.  14,  the  run-time  for  the  full 
model  is  about  seven  times  more  than  the  reduced  model. 

Note  that  the  reduced  model  cannot  be  used  to  calculate  the 
concentration  distribution  inside  the  particles.  If  such  information 
is  valuable,  one  could  use  the  LTI  model  for  most  of  the  particles 
except  for  a  few  selected  ones.  For  those  a  few  selected  particles, 
full  particle  diffusion  equations  are  solved  and  thus  concentration 
distribution  is  available  for  those  a  few  particles.  This  technique 
does  not  apply  to  non-spherical  particles  though  because  it  is  too 


Fig.  11.  Concentrations  of  lithium  ion  in  the  liquid  phase  at  the  four  specified  inter¬ 
faces  obtained  from  the  full  model  (denoted  by  solid  lines)  and  those  from  the  reduced 
model  (denoted  by  dotted  lines)  at  1  C  discharge  rate. 


Fig.  13.  Surface  concentrations  obtained  from  the  full  model  (denoted  by  solid  lines) 
and  those  from  the  reduced  model  (denoted  by  dotted  lines)  at  1  C  discharge  rate  for 
a  few  selected  particles. 
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Fig.  14.  Cell-cycling  simulation  results  based  on  the  full  model  (denoted  by  solid  lines) 
compared  to  results  of  the  reduced  model  (denoted  by  dotted  lines).  Cell-cycling 
protocol:  1  C  discharge  to  3.6  V,  followed  by  1  C  charge  to  4.1  V,  and  then  charge 
the  cell  holding  the  cell  potential  at  4.1  V  until  the  current  decreases  to  0.5  A. 


computationally  expensive  and  too  tedious  to  solve  for  full  3d 
diffusion  equations  even  just  for  one  single  non-spherical  particle 
in  the  physics-based  cell  model.  For  such  non-spherical  particles, 
the  LTI  approach  could  provide  some  additional  information 
without  too  much  additional  computational  cost.  In  the  LTI  models 
generated  so  far,  we  are  only  interested  in  averaged  surface 
concentration.  We  could  have,  though,  added  additional  outputs 
when  creating  the  LTI  model.  For  instance,  we  could  have  used 
concentration  at  a  few  user  specified  locations  inside  the  particle  as 
outputs.  This  would  make  the  LTI  model  have  multiple  outputs 
rather  than  just  one  output.  With  such  an  LTI  model  integrated  into 
the  physics-based  model,  concentration  at  those  selected  locations 
can  be  calculated.  One  of  the  nice  features  about  VF  is  that  adding 
additional  outputs  does  not  increase  the  size  of  the  model  linearly. 
The  VF  method  actually  obtained  its  name  because  of  this  feature. 


decreased  specific  interfacial  area  reflects  the  increased  or 
decreased  surface  area  due  to  different  particle  shapes.  Apart  from 
the  shapes  and  the  specific  interfacial  area,  the  rest  of  the  three 
models  are  identical.  So,  any  difference  in  the  cell  performance  is 
due  to  the  difference  in  particle  shapes.  All  three  models  use  the  LTI 
—  approach  for  the  solid  diffusion  problem. 

1  In  Figs  15  and  16,  discharge  curves  are  compared  for  a  few 

2  different  discharge  rates.  Fig.  15  shows  that  there  is  minimum 
1  difference  using  these  different  particles  when  the  discharge  rates 

are  low  (<1  C  rate).  Fig.  16  shows  that  the  cell  using  the  elliptical 
particles  delivers  more  capacity  at  these  higher  discharge  rates. 
This  is  because  the  surface  area  of  the  elliptical  particle  is  larger  and 
thus  concentration  gradient  inside  the  elliptical  particle  is  smaller 
as  discussed  in  Section  2,  making  it  less  prone  to  saturation  or 
depletion.  On  the  other  hand,  the  cell  using  the  porous  particle  has 
opposite  effect  on  capacity  due  to  its  less  surface  area.  Note  that  the 
porous  particle  here  uses  the  same  spherical  particles  only  with 
overlapping,  so  the  results  showed  the  negative  impact  of  over¬ 
lapping  on  cell  performances. 


5.  Conclusion 

An  LTI  model  is  proposed  for  the  solid-phase  diffusion  problem 
in  particles  used  by  physics-based  lithium-ion  cell  models.  The 
proposed  model  has  the  following  advantages  compared  with 
commonly  used  methods: 

1.  It  can  be  used  for  non-spherical  particles  or  even  porous 
structures.  Loss  of  accuracy  is  minimum  because  the  LTI  model 
can  give  very  similar  results  as  one  would  obtain  by  solving 
a  3d  diffusion  equation. 

2.  The  size  of  the  LTI  model  is  very  small.  For  the  test  cases 
simulated,  the  LTI  model  needs  to  solve  only  7  equations  for 
each  particle,  while  the  full  model  may  need  up  to  200  equa¬ 
tions  for  each  particle  if  spatial  discretization  is  uniform.  So, 
even  for  spherical  particles,  the  proposed  method  can  reduce 
the  total  size  of  the  problem  significantly.  Thus,  the  reduced 
model  runs  much  faster,  a  factor  of  seven  observed  for  the  cell¬ 
cycling  simulation. 


4.  Impact  of  particle  shapes  on  cell  performance 

In  order  to  investigate  the  impact  of  particle  shapes  on  the  cell 
behavior,  three  models  are  created,  which  use  spherical  particles, 
elliptical  particles,  and  porous  particles,  respectively.  The  proper¬ 
ties  of  these  particles  are  listed  in  Table  1.  The  spherical  particles 
have  the  sizes  from  Ref.  [4],  And  the  corresponding  elliptical 
particles  have  the  same  volume.  The  porous  particles  use  the  same 
spherical  particles  with  overlapping.  For  general  porous  structures, 
a  volume  for  a  porous  “particle”  is  not  defined.  So,  the  volume  is 
Table  1  for  the  porous  particle  is  scaled  to  be  per  volume  sense,  and 
its  surface  area  is  also  scaled  to  reflect  that.  The  increased  or 


Particle  properties.  Values  for  the  porous  particle  are  scaled  to  have  the  same 
volume. 


Sphere,  negative 
Sphere,  positive 
Ellipse,  negative 
Ellipse,  positive 
Porous,  negative 
Porous,  positive 


8.1816e— 15 
2.14464e— 15 
8.1816e— 15 
2.14464e— 15 
8.1816e— 15 
2.14464e— 15 


Surface  area 
(m2) _ 

1.96352e—9 
8.0424e— 10 
2.10352e-9 
8.5784e-10 
1.73694e— 9 
7.1143e— 10 


1.0713  xa5,„ 
1.0666X  aSiP 
0.8846X  as  n 
0.8846X  as,p 


The  proposed  LTI  method  greatly  extends  the  capability  of  the 
physics-based  models  by  allowing  for  non-spherical  particles.  And 
the  method  also  significantly  increases  the  efficiency  of  such 
models.  Lastly,  all  are  achieved  with  negligible  loss  of  accuracy. 


0  2,000  4,000  6,000  8,000 

Time  (sec) 


Fig.  15.  Cell  potential  vs.  time  at  different  low  discharge  rates  of  0.5  C  and  1  C  using  the 
three  different  particle  shapes.  No  apparent  differences  due  to  different  shapes  show 
up  at  these  low  discharge  rates. 
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Fig.  16.  Cell  potential  vs.  time  at  different  high  discharge  rates  of  3  C  and  5  C  using 
three  different  particle  shapes. 


List  of  symbols 


as  Specific  interfacial  area,  m  1 

as> n  Specific  interfacial  area  for  negative  electrode,  m  1 

aSrP  Specific  interfacial  area  for  positive  electrode,  m”1 

c  concentration,  mol  m-3 

cs  surface  concentration,  mol  m  3 

D  diffusivity  of  lithium  in  the  solid  particles,  m2  s  1 

j  wall  flux  of  lithium  ions,  mol  m-2  s-1 

L„  thickness  of  negative  electrode,  m 
Lp  thickness  of  positive  electrode,  m 

Ls  thickness  of  separator,  m 

R  radius  of  a  spherical  particle,  m 

t  dimensionless  time 

Appendix  A.  LTI  model  identification  process  for  the  diffusion 
problem  using  vector  fitting 

In  using  the  VF  method  for  state  space  model  identification, 
instead  of  fitting  the  step  response  of  the  state  space  model  to  that 
of  the  CFD  model  directly,  the  transfer  function  of  the  state  space 
model  is  identified  with  the  Fourier  transform  of  the  impulse 
response  of  the  system,  which  is  calculated  from  the  step  response. 

The  transfer  function  of  the  state  space  model  in  Eqn  (4)  can  be 
shown  to  be: 

m  -  Es  qQ.  (A.D 

Transfer  function  is  also  the  Laplace  transform  of  the  impulse 
response  of  the  system.  For  stable  systems,  s  =  ja>  is  in  the  region  of  j= 

convergence  of  the  Laplace  transform  and  thus  fljio)  becomes  the  | 

Fourier  transform  of  the  impulse  response  of  the  state  space  model. 

One  could  then  curve-fit  fijio)  with  the  sampled  Fourier  transform  | 

of  the  impulse  response,  calculated  from  the  step  response,  to  1 

obtain  the  poles  and  residuals  of  Eqn.  (A.1 ).  Such  a  fitting  process  g 

can  be  viewed  as  using  rational  basis  functions  1  /(s  —  a,)  for  curve-  <3 

fitting.  If  dj  were  known,  this  would  become  a  linear  least-squares 
problem  to  determine  q.  Since  a,  are  not  known,  an  iterative 
scheme  is  needed  to  update  a,  at  each  iteration.  Once  the  poles  and 
residuals  are  obtained,  it  is  a  simple  matter  to  obtain  the  corre¬ 
sponding  state  space  model  of  Eqn  (4).  The  sampled  Fourier 
transform  can  be  obtained  by  sampling  the  impulse  response  of  the 


diffusion  problem  followed  by  FFT.  Such  a  fitting  approach  in  the 
frequency-domain  is  generally  referred  to  as  rational  fitting.  VF  is 
a  robust  algorithm  of  rational  fitting  [17],  The  steps  to  identify  the 
state  space  model  using  VF  are  summarized  as  follows: 

Step  1  Obtain  the  step  response  of  the  diffusion  problem  from 
a  CFD  model. 

Step  2  Calculate  the  impulse  response.  The  time  derivative  of  the 
step  response  is  the  impulse  response. 

Step  3  Sample  the  impulse  response  curve. 

Step  4  Perform  FFT  of  the  sampled  impulse  response. 

Step  5  Take  the  low  frequency  portion  of  the  FFT  and  perform 
proper  scaling  in  both  coordinates.  This  gives  the  sampled 
Fourier  transform  of  the  impulse  response. 

Step  6  Perform  vector  fitting  to  obtain  the  poles  and  residuals  of 
the  transfer  function  of  the  state  space  model. 

Step  7  Construct  the  state  space  model  from  the  known  transfer 
function. 

Step  8  Solve  the  state  space  model  in  time-domain. 

Some  variations  of  the  procedure  exist.  First  of  all,  one  could  start 
with  the  impulse  response  if  available.  For  the  diffusion  problem,  it 
is  easier  for  a  CFD  solver  to  calculate  the  step  response  than  the 
impulse  response.  Secondly,  if  the  solver  has  the  frequency  sweep 
capability,  the  sampled  Fourier  transform  can  be  obtained  directly 
and  this  could  then  be  the  starting  point.  The  above  procedure 
assumes  that  the  Fourier  transform  of  the  impulse  response  exists. 
For  the  solid-phase  diffusion  problem,  the  impulse  response  has  no 
Fourier  transfer.  This  is  because  the  surface  concentration 
approaches  to  a  positive  steady  state  value  as  time  goes  to  infinity 
under  an  impulse  flux  boundary  condition.  Such  a  behavior  is  simply 
a  consequence  of  mass  conservation.  In  order  to  apply  the  above 
procedure,  a  small  modification  is  necessary,  which  is  to  apply  the 
above  procedure  after  the  steady  state  value  is  subtracted  from  its 
impulse  response.  And  then  an  integrator  is  added  in  the  final 
transfer  function  to  account  for  the  steady  state  value.  Fig.  Al  shows 
the  numerically  calculated  impulse  response  (not  a  unit  impulse 
response)  for  the  particle  used  in  the  negative  electrode.  The  cor¬ 
responding  step  response  has  a  step  input  of  le— 6  mol  m  2  s-1  into 
the  sphere  rather  than  a  unit  step  input.  Since  the  system  is  linear, 
a  proper  scaling  will  give  the  unit  impulse  response.  Fig.  A2  shows 
the  numerically  calculated  frequency  response  by  performing  FFT  of 
the  impulse  response  followed  by  proper  scaling.  Only  first  1200 
points  from  FFT  is  used  for  fitting  and  displayed  in  Fig.  A2.  Experi¬ 
ence  suggests  that  accurate  results  in  time-domain  can  be  obtained 
when  a  cut-off  frequency  is  chosen  such  that  the  magnitude  of  the 
frequency  response  drops  by  1.5-2  orders  of  magnitude. 


Time  (s) 


Fig.  Al.  Impulse  response  for  the  spherical  particle  for  the  negative  electrode. 
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Fig.  A2.  Numerically  calculated  frequency  response  from  FFT  of  the  impulse  response. 


Appendix  B.  Vector  fitting  results  for  spherical  particles  using 
analytical  solution 

For  spherical  particles,  the  LTI  model  could  be  identified  using 
analytical  step  response  solution  rather  than  CFD  solution.  Such  an 
analytical  solution  in  non-dimensional  form  is  shown  in  Ref.  [8]  to 
be  an  infinite  series, 

C,(T)  =3T+£jQn(T)  (B.l) 

where 

Qn(t)  =4  [l-exp(-^t)]  (B.2) 


A„  -  tan(A„)  =  0  n  =  1,2,... 


(B.3) 


(B.4) 


The  first  500  terms  of  the  series  is  used  here  to  represent  the 
analytical  solution.  Fig.  B1  below  shows  the  step  responses  from  the 
analytical  solution  and  the  LTI  model.  It  can  be  seen  that  very 
accurate  results  are  obtained  using  the  LTI  model.  Below  is  the  state 
space  model  from  the  VF  method  used  to  generate  the  two  plots  in 
Fig.  Bl.  Such  a  state  space  model  can  be  integrated  into  the  physics- 
based  cell  model  for  the  diffusion  problem  using  spherical  particles 
after  scaling  back  to  dimensional  form.  Note  that  matrix  A  in  the 
state  space  model  takes  diagonal  form,  which  makes  the  7 


equations  in  the  state  space  model  decoupled.  Such  a  simple  form 
makes  the  already  small  model  very  efficient  to  solve. 

Initial  conditions:  x0  =  0.0;  xi  =  0.0;  x2  =  0.0;  x3  =  0.0;  x4  =  0.0; 
X5  =  0.0;  x6  =  0.0; 

c_surf  =  366.4550  *  x0  +  38.89143  *  X]  +  16.43058 
*  x2  +  7.834168  *  x3  +  3.825026  *  x4  +  2.198156  *  x5  +  3.0  *  x6; 
dxo/dt  =  -6.084373e  +  004  *  xo  +  fluxj; 
dxi/dt  =  — 5.807664e  +  003  *  Xi  +  fluxj; 
dx2/dt  =  -1.307607e  +  003  *  x2  +  fluxj; 
dx3/dt  =  -3.227365e  +  002  *  x3  +  fluxj; 
dx4/dt  —  -8.353250e  +  001  *  x4  +  fluxj; 
dx5/dt  =  -2.089056e  +  001  *x5  +  fluxj; 
dx6/dt  =  fluxj; 


Fig.  Bl.  Step  responses  in  non-dimensional  form  from  the  analytical  solution  and  the 
LTI  model  for  a  spherical  particle,  a)  Regular  scale  on  i.  b)  Log  scale  on  x. 
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